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MATHEMATICAL  MODELS  OF  NONLINEAR 
GALVANIC  POLARIZATION 


INTRODUCTION 


Serious  corrosion  problems  have  affected  the  hardware  of  the  commercial 
and  military  communities  for  many  years.  It  has  been  estimated  that  a  signif¬ 
icant  percentage  of  the  U.S.  gross  national  product  is  lost  annually  because 
of  ineffective  cathodic  protection.  Corrosion  problems  in  the  commercial  sec¬ 
tor  are  found  in  the  automotive  and  oil  industries;  corrosion  problems  in  the 
military  sector  are  found  in  the  inadequate  protection  of  hulls,  propellers, 
engines,  torpedoes,  missiles,  etc. 

This  study  will  develop  a  numerical  model  that  describes  the  gradual 
corrosion  of  a  surface  from  a  planar  to  a  pitted  contour.  The  model  will 
provide  corrosion  information  such  as  that  found  in  figure  1,  where  pitting 
is  shown  as  a  function  of  time  and  geometry.  The  model  relies  for  its  pre¬ 
dictions  on  established  scientific  laws  (conservation  of  charge  and  Faraday's) 
and  a  fixed  set  of  physical  parameters.  The  parameters  are  chosen  to  provide 
information  on  potential  and  current  distributions  within  the  electrolyte  and 
on  the  anodic  and  cathodic  surfaces.  Sufficient  measurements  are  taken  under 
controlled  conditions  to  ensure  that  the  chosen  electrochemical  parameters 


SIMPLE  ANOO I C/ CATHODIC  GEOMETRY 
ELECTROLYTIC  CONDUCTIVITY  CONSTANT 


o,(t  >  0) 


A  -  ANODIC  SURFACE 
C  -  CATHODIC  SURFACE 
E  -  ELECTROLYTE 
a  -  CONDUCTIVITY 
t  -  TIME 


COMPLEX  ANODIC/CATHODIC  GEOMETRY^ 
UPDATED  ELECTROLYTIC  CONDUCTIVITY 


Figure  la.  Virgin  State  at  t  =  0  Figure  lb.  Evolved  State  at  t  >  0 


Figure  1.  Corrosion  Pitting  as  a  Function  of  Time  and  Geometry 
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result  in  predictions  that  compare  to  reality.  Periodic  potentiostatic  polar¬ 
ization  and  impedance  measurements  are  required  to  update  model  parameters  to 
reflect  accurately  the  state  of  the  pitting  process. 

This  work  integrates  the  geometric  and  electrolytic  evolution  with  time 
into  numerical  models  (finite-element  method)  capable  of  predicting  localized 
corrosion  kinetics.  The  numerical  models  use  empirical  information  to  initiate 
the  predictive  process.  This  report  focuses  on  the  mathematical  treatment  of 
a  general  nonlinear  polarization  layer  consistent  with  the  appropriate  conser¬ 
vation  laws. 


BACKGROUND 


The  characterization  of  localized  corrosion  phenomena,  such  as  pitting 
and  crevice  corrosion,  has  been  studied  for  years.1-3  Mathematical  prediction 
of  electrolytic  corrosion  behavior  for  a  physical  situation  had  its  genesis  at 
least  four  decades  ago.  Traditional  mathematical  techniques  were  reviewed  for 
predicting  the  potential  fields  of  anodes  and  cathodes  in  the  presence  of  geo¬ 
metric  effects,  polarization  curves,  surface  roughness,  etc. 4  Analytical 
attempts  to  calculate  local  current  distributions  in  the  anodic/cathodic 
neighborhood  (surfaces)  demonstrated  the  intractability  of  exact  mathematical 
solutions.5 

In  recent  years,  attempts  have  been  made  to  model  localized  corrosion 
through  use  of  detailed  models  of  the  electrochemistry  in  localized  regions. 
These  models  were  limited  to  simple  geometries  and  constant  electrolytic  prop¬ 
erties.6-9  The  models  were  further  limited  because  the  analyses  did  not  con¬ 
sider  either  geometric  changes  or  changes  in  composition  and  conductivity  of 
the  electrolyte  during  corrosion. 

Recent  studies  have  applied  the  finite-element  numerical  technique  to 
macroscopic  electrogalvanic  field  predictions.  These  models  were  developed  to 
predict  performance  of  cathodically  protected  structures.19-16  Again,  the 
models  did  not  consider  changes  in  geometry  and  electrolytic  properties.  How¬ 
ever,  they  were  successful  in  predicting  macroscopic  current  distributions  at 
the  various  anodic  and  cathodic  areas . 

This  study  investigates  the  use  of  nonlinear  analysis  within  a  subsystem 
of  the  NASA  Structural  Analysis  (COSMIC/NASTRAN)  program  (Heat  Transfer,  Rigid 
Format  3)  to  address  complex  electrode  boundary  conditions  and  electrolytic 
interactions  for  the  beaker  model.  A  parallel  effort,  studying  the  potential 
strength  of  a  more  general  nonlinear  finite-element  program  known  as  MARC,  is 
also  in  progress.16  Both  finite-element  programs  are  maintained  on  a  VAX-11/ 
780  computer  at  the  Naval  Underwater  Systems  Center  (NUSC) . 

A  number  of  papers  demonstrate  that  localized  changes  in  the  electrolyte 
chemistry  occur  with  time  and  geometry  during  corrosion.6  It  is,  therefore, 
important  to  include  these  phenomena  in  any  model  dealing  with  localized  cor¬ 
rosion  kinetics. 
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APPROACH 


The  approach  to  this  research  effort  is  two  phased.  First,  numerical 
models  are  developed  to  describe  the  gradual  corrosion  of  a  surface  from  a 
planar  to  a  pitted  contour.  The  model  provides  corrosion  information  such  as 
that  found  in  figure  1,  where  pitting  is  shown  as  a  function  of  time  and  geom¬ 
etry.  The  model  relies  for  its  predictions  on  established  scientific  laws 
(conservation  of  charge  and  Faraday's)  and  a  fixed  set  of  physical  parameters. 
The  parameters  are  chosen  to  provide  information  on  potential  and  current  dis¬ 
tributions  within  the  electrolyte  and  on  the  anodic  and  cathodic  surfaces. 

In  conjunction  with  this  effort,  sufficient  measurements  are  taken  to 
ensure  the  chosen  electrochemical  parameters  result  in  predictions  that  will 
compare  with  reality.  These  measurements  are  made  under  controlled  laboratory 
conditions  by  Dr.  C.  R.  Crowe,  of  the  Naval  Surface  Weapons  Center  (NSWC) . 
Periodic  potentiostatic  polarization  and  impedance  measurements  are  also 
required  to  update  model  parameters  to  reflect  accurately  the  state  of  the 
pitting  process. 

This  report  emphasizes  the  evaluation  of  existing  finite  elements  for 
nonlinear  applications.  A  number  of  finite  elements  residing  within  COSMIC/ 
NASTRAN  computer  code  are  tested  for  mesh  sizing,  accuracy  of  nonlinear  (table 
input)  membrane  response,  and  various  anodic  and  cathodic  boundary  conditions. 


PHYSICAL  CONSIDERATIONS 


The  introduction  of  mixed-potential  theory  appeared  in  literature  some 
40  years  ago.1?  The  theory  states  that  any  electrochemical  reaction  can  be 
divided  into  two  or  more  oxidation  and  reduction  reactions.  It  negates  the 
existence  of  net-charge  accumulation  during  an  electrochemical  reaction.  The 
significance  of  these  statements  is  that  electrically  coupled  metals  in  an 
electrolyte  behave  in  a  totally  different  manner  than  the  same  metals  in  elec¬ 
trical  isolation. 

The  oxidation-reduction  curves  for  each  metal  are  determined  in  the  labo¬ 
ratory  for  the  steady-state  potential  of  the  metal  in  an  electrolyte.  Figure 
2  demonstrates  this  case  for  zinc,  copper,  and  iron  in  salt  water  (3  percent 
NaCl)  at  30°C  when  measured  in  reference  to  a  standard  electrode,  usually  a 
saturated  calomel  electrode  (SCE)  or  a  silver/silver-chloride  (Ag/AgCl)  elec¬ 
trode.  These  polarization  curves  measure  the  nonlinear  transition  between  two 
potential  states  and  associated  current  density  states .  These  curves  charac¬ 
terize  the  various  electrochemical  reactions  for  an  electrode  independent  of 
exterior  global,  or  large-scale,  conditions.  This  fundamental  assumption 
allows  the  analyst  to  apply  measured  laboratory  uncoupled  electrochemical 
information  to  define  the  electrical  sources  and  sinks  for  a  given  boundary- 
value  problem. 

A  fundamental  assumption  in  this  finite-element  analysis  is  that  the 
base-metal  substrate  has  a  uniform  electrical  potential. The  basis  of  this 
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assumption  is  that  a  base-metal  substrate  has  high  conductivity  compared  to 
that  of  the  electrolyte  and  can,  therefore,  support  large  current  values  for 
small  potential  differences. 

Although  the  metal  substrate  potential  is  assumed  to  be  uniform,  a  poten¬ 
tial  difference  between  the  base-metal  substrate  and  the  electrolyte  of  the 
metal  may  vary  considerably.  (This  difference  is  due  to  paint  or  polarization 
layers  and  can  be  measured  by  a  reference  cell.)  The  potential  difference 
measured  at  a  zinc  surface,  for  instance,  is  certainly  different  from  that 
measured  at  exposed  iron  or  copper  areas. 


THEORETICAL  CONSIDERATIONS 


A  general  finite-element  modeling  (FEM)  .procedure  for  calculating  ele 
galvanic  field  responses  due  to  multiple  anodic/cathodic  interactions  has  a 
developed  for  macroscopic  galvanic  current-density  assessment. 14-15  The  £ 
ic/cathodic  interactions  in  the  conductive  electrolyte  are  predicted  by  th 
application  of  classical  dc  electric-field  theory  for  conductive  continuums  in 
conjunction  with  widely  accepted  laboratory  oxidation/reduction  responses  for 
the  electrodes.  The  electrogalvanic  fields  in  the  electrolyte  are  calculated 
with  the  scalar  Poisson  equation,  whereby  traditional  boundary  conditions  are 
prescribed  in  the  farfield  of  the  electrolyte.  In  the  nearfield  of  the  anodes, 
cathodes,  and  the  painted  metallic  substrate,  complex  boundary  conditions  are 
enforced  based  on  empirical  polarization  curves  and  paint -impedance  values. 

The  ionic  current  in  the  electrolyte,  leaving  the  anode  and  arriving  at  the 
cathodes,  is  mathematically  constrained  to  sum  to  zero  over  the  metallic  sur¬ 
face  (spatial  Kirchhoff's  law). 

Generally,  electric-current  flow  in  a  conducting  medium  is  governed  by 
the  law  of  conservation  of  charge.  In  differential  equation  form,  this  law 
can  be  stated  as 


(1) 


The  constitutive  relationship  (Ohm's  law)  between  current  density  and  the 
electric-field  intensity,  in  terms  of  electrical  conductance,  is 

j  =  oE  ,  (2) 


where,  by  definition,  the  electric-field  intensity  is 

E  *  -V$  .  (3) 


Substitution  of  equations  (3)  and  (2)  into  equation  (1)  yields 


V  • 


3p_ 
3t  ' 


(*0 


This  basic  differential  expression  governing  the  conductive-current  flow 
in  an  electrical  continuum  is  analogous  to  the  general  heat-conduction  equa¬ 
tion.  Specifically,  the  heat-conduction  equation  can  be  stated  as 
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V  •  kVT  +  q  =  pmCh(dT/3t)  .  (5) 

Table  1  illustrates  the  specific  analogous  independent  and  dependent  variables 
relating  equations  (4)  and  (5) .  The  above  analogy  provides  the  analyst  with  a 
mathematical  tool  (FEM)  that  exists  on  numerous  finite-element  computer  codes. 
Many  linear  applications  have  already  been  demonstrated  using  equation  (5) . 

A  need  exists  to  incorporate  mathematically  the  complex  nonlinear  polari¬ 
zations  at  the  various  electrode-electrolyte  interfaces  into  the  field  equa¬ 
tion,  equation  (1);  the  constitutive  equation,  equation  (2);  and  the 
conservation  equation,  equation  (4). 


MODEL  EVOLUTION  AND  DEVELOPMENT 

Various  FEM's  were  generated  to  check  out  (1)  different  elements,  (2) 
nonlinear  responsiveness  by  means  of  table  input  routines,  and  (3)  boundary 
conditions  (single-point  constraint  (SPC) ,  multipoint  constraint  (MPC) , 
enforced  potentials,  and  enforced  currents). 

Figure  3  shows  the  FEM  of  a  steady-state  and  two-dimensional  electrostatic 
problem.  The  model  consists  of  two  four-noded  quadrilateral  elements  with  one 
degree  of  freedom  per  node.  The  electrical  conductivity  of  each  element  is 
assumed  to  be  constant.  The  horizontal  sides  of  the  model  are  assumed  to  be 
insulated.  The  left  nodes  are  held  at  a  constant  electrical  potential.  The 
right  nodes  are  connected  to  another  fixed-potential  state  ’  v  one  of  the 
boundary  elements  under  investigation,  thereby  creating  a  potential  difference. 

Three  different  COSMIC/NASTRAN  elements  were  used  to  model  the  solid-fluid 
interface:  (1)  CELAS2,  (2)  CQDMEM,  and  (3)  CHBDY.  The  CELAS  elements,  spring 
elements,  performed  satisfactorily.  However,  the  element  is  better  suited  for 
structural  than  electrical  applications.  The  electrical  conductivity  must  be 


Table  1.  Comparison  of  Variables 


Thermal 

Variable 

Nomenclature 

Electrogalvanic 

Variable 

K 

Thermal,  electrical  conductivity 

a 

T 

Temperature,  electric  potential 

-VT 

Negative  gradient 

-V<(>  =  E 

-KVT 

Flux,  current  density 

-oV<J>  =  J 

• 

-q 

Negative  time-dependent  source 

(3p/3t) 

Q 

Total  heat,  current  flow 

I 
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Figure  3.  Finite-Element  Model  (FEM)  Used  for 
Boundary- Element  Evaluation 


multiplied  by  a  length  or  area  associated  with  the  boundary  to  obtain  current, 
which  would  be  tedious  for  large,  irregular  meshes. 

It  was  observed  in  an  earlier  computer  run  that  the  CQDMEM  element  readily 
accepted  and  accurately  predicted  the  current  density  from  the  nonlinear 
TABLEM1  card.  At  issue  is  the  correct  geometric  and  impedance  characterization 
of  the  polarization  layer  between  the  metal  surface  and  the  electrolyte.  This 
very  thin  layer  appears  to  control  significantly  the  amount  of  current  density 
and  the  potential  states  (voltage)  at  which  the  anodic  and  cathodic  surfaces 
are  operating. 

It  was,  therefore,  believed  that  the  boundary  element  could  be  modeled  by 
use  of  a  CQDMEM  element  with  extremely  small  thickness  (approximately  10~9 
units) .  However,  if  this  very  thin  quadrilateral  element  was  placed  next  to 
one  of  unit  thickness,  it  would  create  numerical  problems  when  the  conduc¬ 
tivity  matrix  was  formulated.  Therefore,  the  model  was  given  a  logarithmic 
thickness  distribution  (figure  4) .  This  scheme  permitted  use  of  the  CQDMEM 
as  a  boundary  element,  but  it  was  limited  to  a  thickness  of  10"6  units  due  to 
numerical  difficulties.  (The  units  in  these  examples  were  in  inches.)  This 
element  was  not  considered  the  best  choice  because  of  the  many  additional 
elements  and  grid  points  needed  to  model  the  beaker-pitting  situation. 

The  COSMIC  version  of  NASTRAN  possesses  a  boundary  element  called  CHBDY, 
explicitly  intended  for  use  in  modeling  boundary  effects.  It  was  more  attrac¬ 
tive  than  the  other  elements  for  this  application  because  it  did  not  require 
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Figure  4.  Logarithmic  Thickness  Distribution 


additional  nodes  nor  any  correction  to  heat-flux  versus  current-density  values 
The  nonlinear  polarization  trial  function  assigned  to  the  boundary  element  is 
shown  in  figure  5.  The  film  coefficient  was  approximated  by  a  combination  of 
step  functions,  which  were  chosen  because  they  are  representative  of  antici¬ 
pated  potential -dependent  convective  values. 

The  original  configuration  and  computer  results  are  shown  in  figure  6. 

It  was  uncertain  from  these  results  whether  the  program  used  the  average  sur¬ 
face  potential  or  the  average  of  the  surface  and  ambient  potentials  to  obtain 
the  voltage-dependent  film  coefficient. 

The  voltages  of  all  four  boundary  nodes  were  increased  on  the  next  run  to 
ensure  the  program  used  all  potential  ranges  of  the  input  table.  As  figure  7 
shows,  when  the  calculated  average  surface  potential  of  0.5  V  confused  the 
program,  it  used  the  average  value  of  5  V  rather  than  choose  between  3  or  7  V. 
This  run  demonstrated  the  dangers  of  discontinuities  in  nonlinear  property 
tables  and  showed  the  program  accessing  the  table  with  average  surface  voltage 
rather  than  average  film  potential. 

Figure  8  shows  the  last  case,  using  CQDMEM  and  CHBDY  with  increased  sur¬ 
face  potentials.  As  expected,  the  property  table  yielded  a  film  coefficient 
of  7,  and  the  element  conservation  of  current  (equation  (1))  was  maintained. 

It  appears  that  a  single-line  CHBDY  element  does  not  incorporate  the  nonlinear 
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values  from  the  TABLEM1  in  a  satisfactory  manner.  It  is  possible  that  the 
area  CHBDY  element  in  an  eventual  beaker-electrode  model  that  uses  an  area 
element  could  be  more  successful.  This  alternative,  of  course,  must  be  deter¬ 
mined  analytically. 

In  each  of  the  FEM's  mentioned,  the  analysis  was  performed  with  CQDMEM, 
CELAS2,  and  CHBDY  elements.  It  was  difficult  to  determine  whether  the  nonlin¬ 
ear  modeling  problems  were  caused  by  poorly  defined  boundary  conditions, 
improper  reading  by  COSMIC/NASTRAN  of  nonlinear  tables  due  to  lack  of  diag¬ 
nostics,  or  by  numerical  problems  due  to  a  difference  algorithm  within  COSMIC/ 
NASTRAN  (appendix  A,  pages  8.4-1  to  8.4-12,  in  the  NASTRAN  Theoretical  Manual, 
which  addresses  nonlinear  forcing  functions  and  medium  matrices). 18  Conver¬ 
gence  can  be  a  problem  if  some  of  the  numerical  parameters  are  not  properly 
chosen.  However,  a  number  of  guidelines  have  been  suggested  that  generally 
provide  stability  and  convergence  in  the  solution  set.*8 

A  new  solution  technique  was  developed  that  relies  on  the  polarization 
behavior  of  the  given  electrode  material.  For  example,  assume  that  an  anodic 
material  behaves  as  function  f^(<j),i)  in  figure  9  and  that  the  cathodic  elec¬ 
trode  behaves  as  fc($,i),  where  $  is  the  voltage  in  volts  and  i  is  the  current 
density  in  amperes/meter2.  Both  functions  are  determined  empirically.  The 
values  <|>q  and  $A  are  the  constant  open-circuit  potentials  (i  =  0  A)  for  each 
electrode  material.  It  will  be  assumed  that  fA($,i)  and  fc(<f>,i)  were  gene¬ 
rated  within  the  same  electrolyte  under  identical  temperatures  and  conductiv¬ 
ities.  Therefore,  if  these  two  materials  are  forced  to  react  with  each  other, 
the  newly  coupled  system  must  operate  at  some  equilibrium  level  designated  by 
line  A-A' ,  as  indicated  in  figure  9.  Line  A-A'  corresponds  to  a  constant  cur¬ 
rent,  in  which  the  current  leaving  the  anode  must  arrive  at  the  cathode  to 
satisfy  continuity,  if  all  other  boundary  surfaces  are  perfectly  insulated. 
Unfortunately,  the  position  of  line  A-A'  is  not  known  a  priori.  If  a  secant - 
slope  relationship  between  the  potentials  and  the  currents  is  employed,  the 
reduction  curve  for  the  cathode  can  be  defined  and  entered  in  table  2  for  cop¬ 
per.  The  curves  shown  in  figure  9  are  general  representations  for  discussion. 
The  actual  curves  for  copper  are  shown  in  figure  2,19  and  the  conductive  param¬ 
eter,  h(<)0,  is  listed  in  tables  2  and  3  for  the  secant-slope  and  tangent-slope 
definitions,  respectively. 

Polarization  curves  represent  a  relative  voltage  potential  as  a  function 
of  current  density  for  immersed  electrode  materials  in  electrolyte.  Such 
curves  indicate  expected  oxidation/reduction  behavior  for  a  given  electrode 
that  is  subject  to  a  fixed  set  of  conditions  of  the  electrolyte  (temperature, 
degree  of  air  saturation,  composition  of  electrolyte,  etc.).  When  two  elec¬ 
trodes  are  placed  in  an  electrolyte,  mixed  potential  theoryl?  dictates  that, 
for  all  the  partial  oxidation  and  reduction  reactions  that  constitute  an  elec¬ 
trochemical  system,  conservation  of  current,  given  in  equation  (1),  must  be 
maintained.  This  implies  that  point  E  and  point  A  will  shift  in  figure  9  to 
accommodate  the  influence  of  complex  electrochemical  surface  geometries.  The 
true  reaction  for  the  cathodic  surface  lies  somewhere  in  the  neighborhood  of 
point  E,  figure  9.  In  the  horizontal  direction,  the  analytical  solution  to 
the  model  is  physically  constrained  within  the  electrolyte  to  conserve  the 
total  current.  Hence,  the. amount  of  current  leaving  the  anode  must  equal  the 


11 


TR  6921 


Table  2.  Model  of  Copper  Potential  Versus  Current  Curve* 


i  .  .  h( 

(A/m2)  - 


*  Where  1  =  h(*  -  $corr)  and  ^corrosion  =  -°-25  v- 
By  secant-slope  model,  l/h(<$>)  =  ( <j>  -  $Corr)/d- 


!»(♦)  „ 
(A/V  per  m2) 
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amount  of  current  arriving  at  the  cathode.  This  horizontal  bound  is  shown  in 
figure  9  to  be  lines  BCD-B'C'D'. 

The  variation  in  the  potential  state  of  the  cathode  about  point  E  is 
indicated  by  a  lower  secant  slope,  C-C',  and  an  upper  secant  slope,  B-B'.  Each 
secant  line  is  referenced  to  a  differential  potential,  as  indicated  by  the 
intersection  of  line  A-A'  and  lines  B-B'  and  C-C'.  The  secant  slopes  B-B'  and 
C-C'  intersect  at  i  =  0  on  the  vertical  axis. 

By  introducing  a  double-membrane  finite-element  approach,  we  place  a 
bound  on  the  cathodic  operating  potential  (point  E,  figure  9)  of  a  value  above 
(line  B-B')  and  a  value  below  point  E.  Each  membrane  possesses  the  identical 
nonlinear  constitutive  information,  as  measured  in  figure  2.  However,  each 
membrane  operates  at  a  different  potential  value,  although  both  have  the  same 
current  magnitude.  The  average  potential  between  the  two  membranes  is  defined 
as  the  operating  potential • value  for  point  E  (cathodic  operating  potential). 

The  current  flow  from  the  anode  to  the  cathode  is  governed  in  this  approach 
by:  (1)  conservation  of  current,  (2)  upper  cathodic  potential  (open-circuit) 
bound-line  <|>c  =  constant  (figure  9),  (3)  lower  anodic  potential  (open-circuit) 
bound-line  <j>A  =  constant  (figure  9),  and  (4)  the  relevant  surface  dimensions 
and  geometries  causing  the  shift  in  points  E  and  A. 


EXAMPLE  PROBLEM 

The  physical  situation  modeled  ultimately  is  a  small  beaker  partially 
filled  with  an  electrolyte  covering  two  circular  continuous  electrodes,  as 
shown  in  figure  10.  The  area  marked  A  is  the  anode  (iron),  C  is  the  cathode 
(copper),  E  is  the  electrolyte,  and  N  is  the  insulative  material.  The  boundary 
conditions  for  the  beaker  model  are  shown  in  figure  11.  All  the  surface  areas, 
represented  by  3<j>/3n  =  0  (no  current  flow)  are  nonconductive  areas.  The  elec¬ 
trode  surfaces  are  defined  by  means  of  the  double  membrane  explained  in  the 
previous  section.  The  analysis  is  performed  by  an  algorithm  within  NASTRAN 
that  requires  an  initial  guess  for  numerical  iteration. 18  The  mathematical 
example  used  to  verify  the  doubled  nonlinear  membrane  approach  is  shown  in 
figure  12.  This  checkout  model  is  composed  of  two  AREA4  CHBDY  membrane  ele¬ 
ments  spaced  10" 8  m  apart  and  bounded  by  two  spring-type  conductors  (electro¬ 
lyte)  on  the  left  and  a  base  cathodic  potential  on  the  right.  Grid  points  13, 
4,  14,  and  1  represent  the  open-circuit  potential  of  iron  at  -0.65  V  (anodic). 
Scalar  point  99,  on  the  right,  represents  the  base-metal  potential  of  -0.2  V. 
The  conductive  springs  are  modeled  with  a  value  of  4  mho/m,  representing  sea 
water.  The  IR  drops  in  the  FEM  of  figure  12  are  calculated  in  appendix  A. 

The  modeling  results  are  shown  in  appendix  B.  The  total  amount  of  current 
from  the  anode  points  equals  0.6246616  A.  The  amount  of  current  passing 
through  the  membranes  (elements  11  and  12)  and  arriving  at  the  base  metal 
(scalar  point  99)  is  0.6246611  A.  The  continuity  of  current  appears  to  be 
satisfied.  A  similar  analysis  was  completed  by  replacement  of  the  two  spring- 
type  conductors  described  above  with  three-dimensional  finite  elements  (figure 
13) .  The  same  consistency  was  observed  in  the  results  documented  in  appendix 
C. 
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Figure  13.  ONR  CHEXA  Model 


This  approach  should  be  effective  whether  it  is  applied  across  the  anodic 
or  cathodic  electrode  surfaces.  Furthermore,  it  should  not  be  limited  to  a 
specific  number  of  electrodes  if  they  are  interacting  simultaneously. 


RESULTS 


The  example  model  used  a  membrane  finite  element  with  a  cross  section  of 
1  m2.  The  model  (figure  12)  was  assembled  in  units  of  meters.  From  this  prob¬ 
lem  analysis,  the  resultant  current  for  the  defined  system  was  0.62466  A,  which 
corresponds  to  a  potential  in  figure  2  of  approximately  -0.43  V  for  point  E. 

The  potential  value  of  the  first  and  the  second  membrane  are  -0.57191  V  and 
-0.28798  V,  as  shown  in  appendix  B  under  potential  vector.  Because  these  two 
values  represent  the  lower  and  upper  values  of  point  E  (figures  9  and  2) ,  their 
average  indicates  the  actual  potential  state  of  point  E.  That  value  is 
-0.429945  V,  from  the  finite-element  analysis;  it  appears  very  reasonable,  with 
the  nonlinear  correspondence  shown  in  figure  2.  When  the  linear  conductors  in 
figure  12  are  replaced  by  two  three-dimensional  solid  elements  (salt-water 
electrolyte),  as  in  figure  13,  the  current  is  reduced  to  0.4108127  A,  and  the 
resulting  average  potential  at  point  E  is  -0.3512274  V.  For  a  required  shift 
(decrease)  in  system  current  requirements,  a  higher  cathodic  voltage  is 
observed.  This  appears  consistent  with  the  polarization  curve  in  figure  2. 
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CONCLUSIONS 

It  appears  that  the  double  membrane  approximation  for  the  polarization 
layer  is  a  satisfactory  representation  for  nonlinear  electrode  surface  be¬ 
havior.  In  the  mathematical  example  illustrated  in  figures  12  and  13,  a 
secant-slope  impedance  definition  was  used  to  calculate  current  flow  near 
the  electrode  surfaces.  Numerical  convergence  is  realized  for  equation  5  by 
means  of  an  integration  algorithm^  strongly  analogous  to  the  Newmark  B  method 
in  structural  dynamics.  A  tangent  definition  (table  3)  of  impedance  was  also 
used  successfully  for  iron  and  copper  electrodes  in  conjunction  with  a  double 
membrane  configuration. 

A  three-dimensional  FEM  for  the  beaker  will  be  generated  using  the 
double-membrane  approach  in  conjunction  with  three-dimensional  elements 
for  the  electrolyte.  The  beaker  model  will  also  contain  the  effects  of 
nonlinear  anodic  and  cathodic  surface  behaviors  using  both  secant  slope 
and  tangent  slope  definition  for  impedance. 
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Appendix  A 


CALCULATION  OF  IR  DROPS  IN  FEM 


{Base  potential 
=  -0.2  V} 


(table  2) 


Single-dimensional  equivalent 
(from  figure  12) 


$  =  -o.20 

Second  membrane 

7  ~  A<(>  =  0.0879805  cr(<f>) 

i  =  (Ao)A$  =  0.624661 


7  « 


First  membrane 


Water  space 
a=4  mho/m 


•  2 


- «  1 

Fe  surface 
(assume  constant 
<J>  =  -0.65  V) 


o 


3  • 


Acp  =  -0.283936 

i  =  (A(j>)o($)  =  (A<)>)  (2 . 2)  =  0.6746611 
Same  as  element (l) below 


A(|.  =  IR  =  (l/o)  (S./A)  (i) 

Calc  A$  =  (1/4)  (1/1) (0.156165)  =  0.03904 
F.E.  A<t>  =  0.03904 


(a)  Point  E  on  polarization  curve  (figure  2)  was  calculated  to  be  the 
average  of  grid  points  3  and  7. 


-0.571917  -  0.287981 
$E  =  - - - 


-0.429949  V 


(b)  In  equation,  A<|>  =  (l/o)  (i/A)  (i) ,  l  =  1  m,  and  A  =  1  m2. 


A-l/A-2 
Reverse  Blank 
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Appendix  B 

SAMPLE  COMPUTER  PRINTOUT  FOR  ONR  CHUB  ELEMENTS 


This  appendix  provides  a  portion  of  a  typical  computer  printout  of  ONR 
CHUB  elements  for  reference. 


B-l 


MON  NNK.  COPPER  CATHODE  IN  KMCER  MODEL 
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Appendix  C 

SAMPLE  COMPUTER  PRINTOUT  FOR  ONR  CHEXA  ELEMENTS 

This  appendix  provides  a  portion  of  a  typical  computer  printout  of  ONR 
CHEXA  elements  for  reference. 
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INITIAL  DISTRIBUTION  LIST 

Addressee  No.  of  Copies 

Brookhaven  National  Laboratory  (H.  S.  Isaacs)  1 

NAVSEA,  05R25  (H.  Vandervelt) ,  05E1  (W.  Strasburg)  2 

Lehigh  University  (H.  Leidheiser)  1 

DTNSRDC  (B.  Bieberich,  H.  P.  Hack,  W.  Andahazy)  3 

NRl,  6314  (E.  McCafferty)  1 

BK  Dynamics  (J.  Fitzgerald)  1 

ONR,  431  (P.  Clarkin)  1 
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